% This code produces Table 1, 2, 3, and 8. It also produces Figure 1
% and 2.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc

addpath(genpath('Parameter'))
addpath(genpath('Utility'))

addpath(genpath('jplv7'))

%% 0. Import data
Data_sum_file = 'Data/Data_Sum.xlsx';
Data_sum_q = readtable(Data_sum_file,'ReadVariableNames',1,'sheet','Quarterly','basic',true);
Data_sum_pred = readtable(Data_sum_file,'ReadVariableNames',1,'sheet','Pred','basic',true);

dates = gen_eomdate(Data_sum_q.yyyy(1) , Data_sum_q.yyyy(end) , 1 , 12);%1871m1 to 2019m12

l_q = 372;

%% 0.1 Full time series of variables
% For Table 1 and 8
exret = Data_sum_pred.exret;
expd = Data_sum_pred.expd;
expr = Data_sum_pred.expr;
pd = Data_sum_pred.pd_2;% Use repurchase adjusted
infl_ma = Data_sum_pred.infl_ma;
l_q_full = size(Data_sum_pred , 1);

% For Table 2 and 3
exretexp = Data_sum_q.exretexp(end - l_q_full + 1 : end);
retexp = Data_sum_q.retexp(end - l_q_full + 1 : end);
Ret_q = Data_sum_q.Ret(end - l_q_full + 1 : end);
Ret_MS = [nan(3 , 1);
         exp(movsum(log(1 + Ret_q) , 4 , 'Endpoints' , 'discard')) - 1];

% For Table 8
retstd = Data_sum_q.retvol(end - l_q_full + 1 : end);
tb3mreal_ea_BC = Data_sum_q.tb3mreal_ea_BC(end - l_q_full + 1 : end);

%% 1. Table 1

%% 1.1 Parameters
nu_q = 0.018;
Nboot = 50000;
const = 1;
seweight = 1;%NW
block_bs = 1;%1 for i.i.d., 2 for cicular, 3 for stationary

l_q_1 = [372 , 296]; %372 for sample since 1927, 296 for sample since 1946
mainpred_1_name = {'expd_1' , 'expr_1'};
control_1_name = {'infl_ma_1' , 'pd_1'};

Var_RowNames_1 = {'coeff. 1','bias-adj coeff. 1','two-sided p 1', ...
                  'coeff. 2','bias-adj coeff. 2','two-sided p 2', ...
                  'coeff. 3','bias-adj coeff. 3','two-sided p 3', ...
                  'Observations' ,...
                  'R2'};

%% 1.2 Regressions              
Results_1 = cell(length(mainpred_1_name) , length(l_q_1));
mainpred_1 = cell(length(l_q_1) , length(mainpred_1_name));
control_1 = cell(length(l_q_1), 1);
LHS_1 = cell(length(l_q_1), 1);
for i = 1 : length(l_q_1)
    exret_1 = exret(end - l_q_1(i) : end);
    expd_1 = expd(end - l_q_1(i) : end);
    expr_1 = expr(end - l_q_1(i) : end);
    pd_1 = pd(end - l_q_1(i) : end)/100;
    infl_ma_1 = infl_ma(end - l_q_1(i) : end);
    
    nwlags_1 = round(1.3 * (l_q_1(i) - 1)^(1/2));
    
    LHS_1{i} = exret_1;
    
    control_temp = [];
    for q = 1 : length(control_1_name)
        control_temp = [control_temp , eval(control_1_name{q})];
    end
    control_1{i} = control_temp;
    
    for j = 1 : length(mainpred_1_name)
        mainpred_1{i , j} = eval(mainpred_1_name{j});
        RHS_1_temp = [mainpred_1{i , j} , control_1{i}];
        K_1_temp = size(RHS_1_temp , 2);%Number of total RHS variables
        
        res_store_temp_1 = [];
        Var_ColNames_temp_1 = cell(K_1_temp , 1);
        for k = 1 : K_1_temp
            res_temp_1 = pred_VARboot(exret_1 , RHS_1_temp(: , 1 : k) , const , nwlags_1 , seweight , Nboot , block_bs);
            res_vec_temp_1 = [reshape([res_temp_1.b(2 : end) ,...
                              res_temp_1.badj(2 : end) ,...
                              res_temp_1.p2null(2 : end);...
                              nan(K_1_temp - k , 3)]' , [] , 1) ;...
                              l_q_1(i);...
                              res_temp_1.R2badj];
                    
            res_store_temp_1 = [res_store_temp_1 , res_vec_temp_1];  
            Var_ColNames_temp_1{k} = strcat(mainpred_1_name{j}(1 : end - 2) , num2str(k));
        end
        Results_1{j , i} = array2table(res_store_temp_1 ,...
                         'RowNames' , Var_RowNames_1 ,...
                         'VariableNames' , Var_ColNames_temp_1);
    end
end

% Output results
Output_RowNames_1 = {'Experienced payout' , '$\quad$ \{bias-adj. coeff.\}' , '$\quad$ [$p$-value]' ,...
                     'Inflation' , '$\quad$ \{bias-adj. coeff.\}' , '$\quad$ [$p$-value]' ,...
                     '$p-d$' , '$\quad$ \{bias-adj. coeff.\}' , '$\quad$ [$p$-value]' ,...
                     'Observations' , '$R^2$',...
                     'Experienced return' , '$\quad$ \{bias-adj. coeff.\}' , '$\quad$ [$p$-value]' ,...
                     'Inflation' , '$\quad$ \{bias-adj. coeff.\}' , '$\quad$ [$p$-value]' ,...
                     '$p-d$' , '$\quad$ \{bias-adj. coeff.\}' , '$\quad$ [$p$-value]' ,...
                     'Observations' , '$R^2$'};


%% 1.3 Output Table 1
Output_Table_1 = [table2array(Results_1{1 , 1}) , table2array(Results_1{1 , 2}(: , 2 : end));...
                  table2array(Results_1{2 , 1}) , table2array(Results_1{2 , 2}(: , 2 : end))];
              
input_1.data = Output_Table_1;
input_1.tableColLabels = {'1927-2019', '1927-2019' , '1927-2019' , '1946-2019' , '1946-2019'};
input_1.tableRowLabels = Output_RowNames_1;
input_1.dataFormatMode = 'row';

input_1.dataFormat = {'%.2f', 9 ,'%d', 1 , '%.3f' , 1 ,...
                      '%.2f', 9 ,'%d', 1 , '%.3f' , 1};               
input_1.dataNanString = ' '; 
input_1.tableBorders = 0;
latex_Table_1 = latexTable(input_1 , 'Results/Table_I.tex');

%% 2. Figure 1

%% 2.1 Future realizations
cum_length = 20; %5-year

exret_lr_data = cell(length(l_q_1) , 1);

%% 2.2 Predicted
exret_lr_fitted = cell(length(l_q_1) , length(mainpred_1_name));
for i = 1 : length(l_q_1)
    LHS_temp = LHS_1{i};
    exret_lr_data{i} = movsum(LHS_temp(2 : end) , cum_length , 'Endpoints' , 'discard');
    for j = 1 : length(mainpred_1_name)
        mainpred_temp = mainpred_1{i , j};
        
        % AR(1) coefficients for the main predictor
        arcoeff_temp(2) = 1 - nu_q;
        arcoeff_temp(1) = mean(mainpred_temp) * (1 - arcoeff_temp(2));
        
        % Predictive coefficients
        pred_temp = pred_VARboot(LHS_1{i} , mainpred_temp , const , nwlags_1 , seweight , Nboot , block_bs);
        
        % n-period ahead forecasts
        mainpred_ts = bsxfun(@times , mainpred_temp(1 : end - 1) , arcoeff_temp(2).^(0 : cum_length - 1)) + repmat(arcoeff_temp(1) * (1 - arcoeff_temp(2).^(0 : cum_length - 1))/(1 - arcoeff_temp(2)) , l_q_1(i) , 1);
        
        % Sum of n-period ahead forecasts
        exret_lr_fitted{i , j} = sum(mainpred_ts * pred_temp.badj(2) + pred_temp.badj(1) , 2);
    end
end

%% 2.3 Output Figure 1
exret_lr_fitted_expd = exret_lr_fitted{1 , 1};

figure
plot(dates(end - 3 * l_q_1(1) + 3 : 3 : end - 3 * cum_length + 3) , exret_lr_data{1} * 100 , 'b' , 'LineWidth' , 2 )
hold on
plot(dates(end - 3 * l_q_1(1) + 3 : 3 : end - 3 * cum_length + 3) , exret_lr_fitted_expd(1 : end - cum_length + 1) * 100 , 'r' , 'LineWidth' , 2 )
ylim([min(exret_lr_data{1}) 1.5 * max(exret_lr_data{1})] * 100)
recessionplot
legend('Actual future 5-year excess returns',...
       'Predicted future 5-year excess returns' , 'FontSize' , 12)
ylabel('Cumulative log excess returns %' , 'FontSize' , 12)
xlabel('Date' , 'FontSize' , 12)
set(gca,'TickLabelInterpreter', 'latex')
saveas(gcf,'Figures/Figure1','epsc')
savetightfigure('Figures/Figure1')
    
%% 4. Table 2

%% 4.1 Panel A regressions
% We regress survey expected excess returns in quarter t on experienced
% growth in quarter t-1 (to make sure information is known to households)
LHS_2A = exretexp(2 : end);%Subjective expected excess returns

lexpd_2A = expd(1 : end - 1);
lexpr_2A = expr(1 : end - 1);
control_2A = Ret_MS(1 : end - 1);%We use lagged one-year return (in percentage) for control

mainpred_2A_name = {'lexpd_2A' , 'lexpr_2A'};

Var_RowNames_2A = {'coeff. 1','s.e. 1', 'p-value 1' ,...
                   'coeff. control','s.e. control', 'p-value control', ...
                   'Const','s.e. Const', 'p-value Const', ...
                   'Observations' ,...
                   'R2adj'};               
            
Results_2A = cell(length(mainpred_2A_name) , 1);
for j = 1 : length(mainpred_2A_name)
    mainpred_2A_temp = eval(mainpred_2A_name{j});
    RHS_2A_temp = [mainpred_2A_temp , control_2A];
    K_2A_temp = size(RHS_2A_temp , 2);%Number of total RHS variables
    
    res_store_temp_2A = [];
    Var_ColNames_temp_2A = cell(K_2A_temp , 1);
    for k = 1 : K_2A_temp
        % Estimate
        [b_temp_2A , se_temp_2A , p_temp_2A , Radj_temp_2A , T_temp_2A] = olsnan_gmm_EWC(LHS_2A , RHS_2A_temp(: , 1 : k) , 1);%Calculate s.d. using EWC
        % Store results
        res_vec_temp_2A = [reshape([b_temp_2A(2 : end) , se_temp_2A(2 : end) , p_temp_2A(2 : end);...
                                    nan(K_2A_temp - k , 3) ;...
                                    b_temp_2A(1) , se_temp_2A(1) , p_temp_2A(1)]' , [] , 1) ;...
                                    T_temp_2A;...
                                    Radj_temp_2A];       
                                  
        res_store_temp_2A = [res_store_temp_2A , res_vec_temp_2A];  
        
        Var_ColNames_temp_2A{k} = strcat(mainpred_2A_name{j}(2 : end - 3) , num2str(k));
    end
    Results_2A{j} = array2table(res_store_temp_2A ,...
                               'RowNames' , Var_RowNames_2A ,...
                               'VariableNames' , Var_ColNames_temp_2A);                      
end
                                                            
%% 4.3 Panel B regressions
LHS_2B = Ret_MS(6 : end) ... % Because control ends in t-1 and realization ends in t+4
       - retexp(2 : end - 4);
   
% We regress survey expected excess returns in quarter t on experienced
% growth in quarter t-1
control_2B = Ret_MS(1 : end - 5);%We use lagged one-year return (in percentage) for control
lexpd_2B = expd(1 : end - 5);
lexpr_2B = expr(1 : end - 5);
   
mainpred_2B_name = {'lexpd_2B' , 'lexpr_2B'};

Var_RowNames_2B = {'coeff. 1','s.e. 1', 'p-value 1'...
                   'coeff. control','s.e. control','p-value control' ...
                   'Const','s.e. Const' , 'p-value Const', ...
                   'Observations' ,...
                   'R2adj'};   
               
Results_2B = cell(length(mainpred_2B_name) , 1);
for j = 1 : length(mainpred_2B_name)
    mainpred_2B_temp = eval(mainpred_2B_name{j});
    RHS_2B_temp = [mainpred_2B_temp , control_2B];
    K_2B_temp = size(RHS_2B_temp , 2);%Number of total RHS variables
    
    res_store_temp_2B = [];%s.d. calculated using NW
    Var_ColNames_temp_2B = cell(K_2B_temp , 1);
    for k = 1 : K_2B_temp
        % Estimate
        [b_temp_2B , se_temp_2B , p_temp_2B , Radj_temp_2B , T_temp_2B] = olsnan_gmm_EWC(LHS_2B , RHS_2B_temp(: , 1 : k) , 1);%Calculate s.d. using EWC       
        % Store results
         res_vec_temp_2B = [reshape([b_temp_2B(2 : end) , se_temp_2B(2 : end) , p_temp_2B(2 : end);...
                                      nan(K_2B_temp - k , 3) ;...
                                      b_temp_2B(1) , se_temp_2B(1) , p_temp_2B(1)]' , [] , 1) ;...
                                      T_temp_2B;...
                                      Radj_temp_2B];                       
                    
        res_store_temp_2B = [res_store_temp_2B , res_vec_temp_2B];  

        Var_ColNames_temp_2B{k} = strcat(mainpred_2B_name{j}(2 : end - 3) , num2str(k));
    end
    Results_2B{j} = array2table(res_store_temp_2B ,...
                               'RowNames' , Var_RowNames_2B ,...
                               'VariableNames' , Var_ColNames_temp_2B);
end
                                         
%% 4.4 Output Table 2   
input_2.tableColLabels = {'(1)', '(2)' , '(3)' , '(4)'};

input_2.tableRowLabels = repmat({'Experienced' , ' ' , ' ', ... 
                                 'Lagged one-year return' , ' ' , ' ',...
                                 'Constant' , ' ' , ' ',...
                                 'Observations' , 'Adj. $R^2$'} ,...
                                 1 , 2);

input_2.data = [table2array([Results_2A{1}(: , 1) , Results_2A{1}(: , 2) ,Results_2A{2}]) ;...
                table2array([Results_2B{1}(: , 1) , Results_2B{1}(: , 2) ,Results_2B{2}])];

input_2.dataFormatMode = 'row';
input_2.dataFormat = {'%.2f', 9 ,'%d', 1 , '%.3f' , 1 ,...
                      '%.2f', 9 ,'%d', 1 , '%.3f' , 1};   
                  
input_2.dataNanString = ' '; 
input_2.tableBorders = 0;

latex_Table_2 = latexTable(input_2 , 'Results/Table_II.tex');

%% 5. Table 3

%% 5.1 Regressions
realltg = Data_sum_q.realltg(end - l_q_full + 1 : end);
realmtg = Data_sum_q.realmtg(end - l_q_full + 1 : end);

LHS_3 = realltg(2 : end);

lexpd_3 = expd(1 : end - 1);
lexpr_3 = expr(1 : end - 1);
control_3 = [realmtg(2 : end) , Ret_MS(1 : end - 1)];%We use lagged one-year return (in percentage) for control

mainpred_3_name = {'lexpd_3' , 'lexpr_3'};

Var_RowNames_3 = {'coeff. 1','s.e. 1','p 1', ...
                  'coeff. control 1','s.e. control 1', 'p control 1', ...
                  'coeff. control 2','s.e. control 2', 'p control 2' ...
                  'Const','s.e. Const', 'p Const 2'...
                  'Observations' ,...
                  'R2adj'};

Results_3 = cell(length(mainpred_3_name) , 1);
for j = 1 : length(mainpred_3_name)
    mainpred_3_temp = eval(mainpred_3_name{j});
    RHS_3_temp = [mainpred_3_temp , control_3];
    K_3_temp = size(RHS_3_temp , 2);%Number of total RHS variables
    
    res_store_temp_3 = [];
    Var_ColNames_temp_3 = cell(K_3_temp , 1);
    for k = 1 : K_3_temp
        % Estimate
        [b_temp_3 , se_temp_3 , p_temp_3 , Radj_temp_3 , T_temp_3] = olsnan_gmm_EWC(LHS_3 , RHS_3_temp(: , 1 : k) , 1);%Calculate s.d. using EWC       
        % Store results
        res_vec_temp_3 = [reshape([b_temp_3(2 : end) , se_temp_3(2 : end) , p_temp_3(2 : end);...
                                   nan(K_3_temp - k , 3) ;...
                                   b_temp_3(1) , se_temp_3(1) , p_temp_3(1)]' , [] , 1) ;...
                                   T_temp_3;...
                                   Radj_temp_3];   
                                 
        res_store_temp_3 = [res_store_temp_3 , res_vec_temp_3];  
        
        Var_ColNames_temp_3{k} = strcat(mainpred_3_name{j}(2 : end - 2) , num2str(k));
    end
    Results_3{j} = array2table(res_store_temp_3 ,...
                               'RowNames' , Var_RowNames_3 ,...
                               'VariableNames' , Var_ColNames_temp_3);
end


%% 5.2 Output Table 3
input_3.tableColLabels = {'(1)', '(2)' , '(3)' , '(4)' , '(5)' ,'(6)'};

input_3.tableRowLabels = {'Experienced' , ' ' , '' , ... 
                          'One-year EPS growth forecast' , ' '  , ' ' ,...
                          'Lagged one-year return' , ' ' , ' ' ,...
                          'Constant' , ' ' , ' ' ,...
                          'Observations' , 'Adj. $R^2$'};

input_3.data = table2array([Results_3{1} , Results_3{2}]);

input_3.dataFormatMode = 'row';
input_3.dataFormat = {'%.2f', 12 ,'%d', 1 , '%.3f' , 1};   
                  
input_3.dataNanString = ' '; 
input_3.tableBorders = 0;

latex_Table_3 = latexTable(input_3 , 'Results/Table_III.tex');

%% 6. Figure 2
realltg_nm = realltg(nonnan_first(realltg) : end);

figure
plot(dates(end - 3 * length(realltg_nm) + 3 : 3 : end) , realltg_nm * 100 , 'b' , 'LineWidth' , 2)
hold on
plot(dates(end - 3 * length(realltg_nm) + 3 : 3 : end) , expd(end - length(realltg_nm) : end - 1) * 100 , 'r' , 'LineWidth' , 2)
ylim([-0.1 , 4])
legend('Long-term earnings growth forecast' , 'Experienced real payout growth' , 'FontSize' , 12)
xlabel('Date' , 'Fontsize' , 12)
ylabel('Growth rates %' , 'Fontsize' , 12)
set(gca,'TickLabelInterpreter', 'latex')
saveas(gcf,'Figures/Figure2','epsc')
savetightfigure('Figures/Figure2')

%% 7. Table 8

%% 7.1 Regressions
l_q_8 = [372 , 296]; %372 for sample since 1927, 296 for sample since 1946

Results_8 = cell(length(l_q_8) , 1);
for i = 1 : length(l_q_8)
    LHS_8 = retstd(end - l_q_8(i) : end).^2;
    expd_8 = expd(end - l_q_8(i) : end);
    pd_8 = pd(end - l_q_8(i) : end)/100;
    tb3mreal_ea_8 = tb3mreal_ea_BC(end - l_q_8(i) : end);
    
    [b_1 , se_1 , p_1 , Radj_1] = olsnan_gmm_EWC(LHS_8(2 : end) , expd_8(1 : end - 1) , 1);
    [b_2 , se_2 , p_2 , Radj_2] = olsnan_gmm_EWC(LHS_8(2 : end) , pd_8(1 : end - 1) , 1);
    [b_3 , se_3 , p_3 , Radj_3] = olsnan_gmm_EWC(LHS_8(2 : end) , [pd_8(1 : end - 1) , tb3mreal_ea_8(1 : end - 1)] , 1);
    
    % Store results                
    Results_8{i} = [b_1(2) , NaN , NaN ;...
                    se_1(2) , NaN , NaN ;...
                    p_1(2) , NaN , NaN ;...
                    NaN , b_2(2) , b_3(2) ;...
                    NaN , se_2(2) , se_3(2) ;...
                    NaN , p_2(2) , p_3(2) ;...
                    NaN , NaN , b_3(3) ;...
                    NaN , NaN , se_3(3) ;...
                    NaN , NaN , p_3(3) ;...
                    b_1(1) , b_2(1) , b_3(1) ;...
                    se_1(1) , se_2(1) , se_3(1) ;...
                    p_1(1) , p_2(1) , p_3(1) ;...
                    l_q_8(i) , l_q_8(i) , l_q_8(i);...
                    Radj_1 , Radj_2 , Radj_3];            
end

%% 7.2 Output Table 8
input_8.tableColLabels = {'(1)', '(2)' , '(3)' , '(4)' , '(5)' ,'(6)'};

input_8.tableRowLabels = {'Experienced real payout growth' , ' ' , ' ' , ... 
                          '$(p - d)/100$' , ' ' , ' ' ,...
                          'Real risk-free rate' , ' ' , ' ' ,...
                          'Constant' , ' ' , ' ' ,...
                          'Observations' , 'Adj. $R^2$'};

input_8.data = [Results_8{1} , Results_8{2}];

input_8.dataFormatMode = 'row';
input_8.dataFormat = {'%.2f', 12 ,'%d', 1 , '%.3f' , 1};   
                  
input_8.dataNanString = ' '; 
input_8.tableBorders = 0;

latex_Table_8 = latexTable(input_8 , 'Results/Table_VIII.tex');
  
